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Abstract 

A novel approach for simulating potential propagation in neuronal branches with high accuracy is developed. 
The method relies on high-order accurate difference schemes using the Summation-By-Parts operators with 
weak boundary and interface conditions applied to the Hodgkin-Huxley equations. This work is the first 
demonstrating high accuracy for that equation. Several boundary conditions are considered including the 
non-standard one accounting for the soma presence, which is characterized by its own partial differential 
equation. Well-posedness for the continuous problem as well as stability of the discrete approximation is 
proved for all the boundary conditions. Gains in terms of CPU times are observed when high-order operators 
are used, demonstrating the advantage of the high-order schemes for simulating potential propagation in large 
neuronal trees. 

Keywords: High-order accuracy; Hodgkin-Huxley; Neuronal networks; Stability; Summation- by-parts; 
Well-posedness 



1. Introduction 

Understanding the integration of synaptic input by a neuron and the propagation of the signal to its 
own output synapses is of high importance in neurosciences. Numerical simulation of such a phenomena 
has become an option since Hodgkin and Huxley developed their model in 1952 [1 j. The Hodgkin-Huxley 
equations are a set of coupled partial and ordinary differential equations. The first one is the cable equation 
that describes the distribution and evolution of the intracellular potential. The other equations are related 
to the evolution of gating variables describing ion channels dynamics inside a neuron, which is typically 
constituted of dendritic branches, a soma, an axon and synapses. Appropriate boundary conditions can be 
associated with branch ends, junctions and the soma. The soma boundary condition is non-standard since 
it consists of a linear relation between time and space derivatives of the solution and the solution itself. 

Subsequently, Rail developed methods for solving the potential propagation in passive neuronal trees 
where branches satisfying the "3/2" law are connected to the soma [2 . Hines then simulated the potential 
propagation in a neuronal tree without soma using second order finite difference schemes [3]. That scheme 
was subsequently used for predictions of potential propagation in dendritic trees [4j |5] . 

In the present work, high-order accurate difference schemes based on Summation-By-Parts (SBP) opera- 
tors with the weak Simultaneous Approximation Term (SAT) procedure 151 flQ| rTT | rT2| fT3 l fT4l fT5 | fTB] 
are applied to the Hodgkin-Huxley equations and their associated boundary conditions. Rigorous conver- 
gence properties are demonstrated, even in the presence of the non-standard soma boundary condition. The 
ability to apply high order schemes for the solution of the Hodgkin Huxley equations is essential since it 
results in a lower computational cost for the large systems that arise when neurons with large dendritic trees 
are considered. It will be demonstrated that the SBP-SAT technique provides a modular approach that 
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is particularly effective for coupling of branches in a dendritic tree. Indeed, the SBP-SAT method offers a 
very effective way to enforce potential continuity and current conservation at the junction between those 
branches. 

This paper is organized as follows: in Section [2j the continuous set of partial differential equations and 
associated boundary conditions are presented and their strong well-posedness is demonstrated using energy 
estimates for two essential cases: (1) the case of an axon connected to the soma and (2) the case of a 
dendritic tree. In Section |3j discretization of each of the aforementioned problems is carried out and the 
associated penalty coefficients are chosen so that semi-discrete energy estimates hold. This results in identical 
estimates as in the continuous case and strong stability of the semi-discrete problem. In Section [4j the order 
of convergence associated with each of the two aforementioned problems is investigated using the method of 
manufactured solutions. Applications of the new proposed approach to large neuronal systems are reported 
in Section [U Conclusions are drawn in Section [U 



2. The continuous problem 

2.1. Equations 

The Hodgkin-Huxley equations [1 are a set of coupled partial and ordinary differential equations ex- 
pressed in terms of the (1) intracellular potential u(x,t) and (2) gating variables m(x,t), h(x,t) and h(x,t) 
describing ion channels dynamics. The computational domain is in this section x G [0,L]. 

The equation for the potential u is based on the cable equation [17j and can be written as 

ut = -j^(a(x) 2 u x ) x - -^-g(m(x,i),h(x,i),n(x,t))u 
CiyX) O m 

+ -^f(m(x,t),h(x,t),n(x,t),x,t), (x,t) e [0,L] x [0, T]. 
In 0, a(x) is the radius of the neuron at the location C m the specific membrane capacitance and 

where Ri denotes the axial resistivity. The conductance g(m, ft, n) of the cable in terms of the gating variables 
is 

g(m, ft, n) = gim 3 h + g 2 n 4 + g 3 , (2) 
where gi > 0, i = 1, 2, 3. The expression for /(ra, ft, n, x, t) is given by 

f(m, ft, n, x, t) = giEim 3 h + g2E 2 n 4 + g 3 E 3 - I(x, i), (3) 

where i = 1, 2, 3 are equilibrium potentials. I(x,t) is an input current at location x that originates either 
from artificial current injection or synaptic input from another neuronal cell. 
The equations describing the evolution of the gating variables are 

m t (x,t) = a m (u(x,t))(l - m(x,t)) - f} m (u(x,i))m(x,i), 

ht(x,i) = a h (u(x,i))(l - h(x,i)) - p h (u(x,t))h(x,i), (4) 
n t (x,t) = a n (u(x,t))(l — n(x,t)) — /3 n (u(x,t))n(x,t). 

where (x,t) G [0,L] x [0,T]. Expressions for a m , a^, a n , /3 m , fa and fa are provided in Appendix A. 

One can prove the following property based on the equations associated with the gating variables. 
Proposition 1. Let x G [0,L]. If < m(x,0) < 1 and if m(x,i) is C° in time, then < m{x,t) < 1 \/t G 
[0,21- 

Proof: see Appendix B. 
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Figure 1 Outer normal vectors to a branch 



The same property also applies for h and n. As a result, the function g can be bounded as follows 

3 

< #3 < g(m,n,h) = gim 3 h + g 2 n 4 + g 3 < (5) 

2 = 1 

The bound Q will be useful when well-posedness of the initial boundary value problems is studied in 
Section 12.61 

Several types of boundary conditions can be associated with the cable equations [17] . 

2.2. Voltage clamp boundary conditions 

The simplest boundary condition associated with Eq. ([I]) is the prescribed potential condition at x& = 
or x b = L: 

u(xb,i) = u (t). (6) 
This corresponds to controlling the potential at the extreme end of the cable. 

2.3. Sealed end boundary conditions 

Another possible boundary condition for Eq. is the sealed end boundary condition at x b £ {0,L}: 

Vu(x h ,t) ■ n(x b ) = 0, (7) 

where n(x b ) is the outer normal vector to [0, L] at the end point x b of the cable, as depicted in Figure [I] 
Vu(x bl i)-n(xb) corresponds, up to a constant factor, to the current that exits the domain [0, L] at x b . Hence, 
enforcing the sealed end boundary conditions is equivalent to prescribing that there is no current exiting the 
neuron at x b . 

Remark. Although the field u only depends on one space variable, in the remainder of the paper, the 
notation \7u(x b ,t) • n(xb) is used in place of —u x (xb,t) when x b is a left boundary (that is n(xb) = —1) and 
in place of u x {x^t) when x^ is a right boundary (n(xb) = 1). 

2.4- Soma boundary conditions 

When a soma is located at an extremity x^ of a cable, a boundary condition describing the current 
conservation in the soma applies. This configuration is graphically depicted in Figure [2] (a). The boundary 
condition is 

ut(xb, t) = -rja(x b ) 2 Vu(x b , t) • n(x b ) - -^-g(m(x b , i), n(x h , t), h(x b , t))u{x h , t) 

! m (8) 

+ y^f( m ( x b, t), n(x b , i), h(x h , t),x b ,t) 

with 

where A soma denotes the soma surface area. 

Remark. Note here the similarity between the boundary equation ([8| and the cable equation ([!]). The only 
difference lies in the spatial derivative term. 
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(a) Cable and soma (b) Junction 

Figure 2 Two configurations for a neuronal network 



2.5. Junction interface conditions 

When multiple cables meet at a common junction, as depicted in Figure [2] (b) , appropriate interface con- 
ditions have to be applied. These conditions correspond to potential continuity and current conservation [4]. 
If N c cables of respective potential wW, i = 1, • • • , N c and cable radius a^(-) have a junction at x = x b: 
then these interface conditions are 

uW(x b ,t)=uW(x b ,t), i,je{l r -,N c } (10) 

for potential continuity and 

N c 

£(a« (x h )) 2 Vu^ (x bj t) • ra« (x b ) = (11) 

i=l 

for current conservation. 



2. 6. Well-posedness 
2.6.1. Cable and soma 

The case of a single cable connected to the soma is first studied. Dropping the source terms, the cable 
equation is written as 

a(x)u t = ii(a(x) 2 u x ) x - °^g(x, t)u, (x, t) e [0, L] x [0, T] (12) 
where g(x,i) = g(m(x,t),n(x,t),h(x,t)) e bmin,5max] and a(x) e [a min ,a max ], with 

3 

#min = 93 > 0, # max = ^9%, «min > 0. (13) 

The boundary conditions associated with the cable equation are, in the case of a single cable connected to 
the soma at x b = L, 

1^(0) = 0, u t (L) = -rja(L) 2 u x (L) - -±-g{L,t)u(L). (14) 
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Defining the energy norm of u at time t as 



Ml 2 



u(x, t) 2 dx. 



By multiplying Eq, (12) with u and integrating over the domain [0,L], one obtains 



;||vM? = -llv^^ll 2 + »(a(L) 2 u x (L)u(L) -a(0) 2 u x (0)u(0)) 



The sealed boundary condition at x = and the soma boundary equation at x = L result in 




1 



||VSti||? + £ (u 2 ) t (L) = -\\^au x \\ 



2" v 2r? 
Defining the weighted norm 



2 WCM) 



Hl2 = ||V5ti|| 2 + ^(i) 2 , 




(15) 



(16) 



(17) 



Eq. (17) becomes, 



(\\u\\t) t + 2\\ y ^au x \\ 2 + 2 
This result can be summarized in the following Proposition. 



0. 



(18) 



Proposition 2. The initial boundary value problem (12)-(14) is well-posed. 



2.6.2. Cables with junction 

In this section, the case of N c cables connected at a junction is considered. For simplification, it is 
assumed that all the other cable extremities are either sealed or clamped. If one of those extremities was 
connected to the soma, a similar approach to the one presented in the previous section can be used. 

(x®,t) G [0,LW] x [0,T] here denotes the potential in the z-th cable, i = l,--- ,N C . To 
simplify the notations, it is assumed here that the junction is located at x^ = 0, all the other cable 
extremities being sealed, except the last one where the potential is held constant at zero. Hence, the PDEs 
of interest are of the same form as ( 12 ) where u : x, a and L are replaced by ? #0) ? a M anc [ £,M respectively. 
The interface conditions are 



M W(0,t)=ti^(0,t), Vi,j = l,-..,JV C , ^(a^(0)) 2 ^ } (0,t)=0, 

i=l 

ug 4) (L«t) = 0, £ = 1,...,JV C -1 



and 



(19) 

(20) 
(21) 

For each cable i G 1, • • • , iV c equation ([16]) holds with appropriate variables. Summing those equations 



,W0 



(L^),t) = 0. 



The energy norm for at time t is defined in ([15]) where the integral bounds are and LW. 



and applying the sealed end and voltage clamp boundary conditions leads to 
2 



1 iV c 2 iV c 2 AT C 



= 1 



= 1 



E 



=i 



(22) 
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Applying the potential continuity condition at the junction, u^\0) = u^(0), i = 2, • • • , iV c , the identity 



^(a«(0)) 2 ^ ) (0) U «(0) = W ( 1 )(0)^(aW(0)) 2 ^ ) (0) 



(23) 



i=l i=l 

and the current conservation at the junction leads to cancelation of the second factor in the right hand side 



of equality (22). Hence 



1 ^ c II / 2 

Ij2 v^ (i) 



i=l 



One can define the weighted norm || • ||* for u = [u^\ • • • , u^ Nc ^] T as 



h\\i = J2 



1=1 



The energy equality then becomes 



(\\u\\l) t + 2\\^mu x \\l+2j^, 



j(0oW 



0. 



This result is summarized in the following Proposition. 



(24) 



(25) 



(26) 



Proposition 3. The initial boundary value problem (12), (19), (20), (21) holding for N c cables connected 



at a junction and with other cable extremities sealed or clamped is well-posed. 



3. The semi-discrete problem 



In this section, Eq. ( |12| ) is discretized on the domain [0, L] using a uniform mesh of N + 1 points. The 
discrete approximation of u(-,i) is 



u(t) = [u (t), • • • , u N (t)] T , Ui(t) « u 



i - 1 



L,t , i = 1, 



,iV. 



Operators of SBP form [18 approximate the derivative of u(-,i) as 

u x (t) = [u x0 (t), • • • , ^iv (t)] T w Diu(t) = P^QuW, 



(27) 



where P is a diagonal symmetric positive definite matrix and Q satisfies Q + Q T = B = diag(— 1, 0, • • • , 0, 1). 
The second space derivative in this paper is approximated as 

u xx (t) * D 2 u(t) = DiDiu(t) = (P- X Q) 2 u(t). (28) 

Compact second derivatives in second order form also exist [TOJ [19] . 

To alleviate notations, the dependency of u and its derivatives on time is omitted in the following. A 
norm based on P is 

||u|| = Vu T Pu. (29) 



The semi-discrete version of ( 12 ) without inclusion of boundary conditions is 



Au, = /xP-^AV) - — AG(t)u 



(30) 



where the diagonal matrix A = diag(ao, • • • , a/v) an d the diagonal matrix G(t) = diag(g (£), • • • > #iv(£))- 
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3.1. Cable with soma 



The boundary conditions ([14]) are integrated in the formulation as penalty terms in the following way 

(31) 



Aut = mP _1 Q(A 2 u x ) - — AG(t)u + troP-^Uxo - 0)e 



a L P 1 ( u tN + V a N u xN + ^-9N{t)u N - ) e L , 



where, for clarity, the data is indicated by 0. 
Premultiplying Eq. (31) by u T P leads to 



u T PAu t = /iu T Q(A 2 u x ) - - r u T PAG(t)u + a u x0 u 

+ 0~LU N [UtN + V a N u xN + ~^-9N{t)u N 

Since P, A and G(t) are diagonal matrices, they commute and lead to 
1 



T^.^r x T G(M„ G(t)A 



/ G(t)A 



u T PAu t = u T v/AP>/Au t = illv'Aull, 2 . 



The SBP properties lead to 



u T Q(A 2 u,) = u T (B - Q t )P- 1 P(A 2 u,) 

= -u t (P- 1 Q) t P(A 2 u.) + u T B(A 2 u.) 
= -u^P(A 2 u x ) - clIuqu x q + a 2 N u N u xN 
= -||Au x || 2 + a 2 N u N u xN - alu u x0 . 



By the choice ctq = /iafj and cr^ = — ^, Eq. (32) becomes 



Aull 



-fi\\Au x 



'G(t)A 



/i / 1 

: ^AT ( WtJV + -^-9N(t)U N 



Similarly as in the continuous case, one can define a weighted norm 



Au 

which leads following energy estimate and proposition 

||u||2 + 2 M ||A U;c || 2 + 2 



'G(t). 



C77 



= 0. 



(32) 



(33) 



(34) 



(35) 



(36) 



(37) 



(38) 



Proposition 4. The SBP-SAT scheme (31) for solving the semi-discrete problem associated with the cable 
equation with soma as a boundary condition is stable [18j |20] . 



mate (18) derived in Section 2.6.1 



Remark. Note that the energy estimate (38) is the semi-discrete analog to the continuous energy esti- 
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3.2. Cables with junction 

Each cable is here discretized on the domain [0, L^], i = 1, • • • , N c using a uniform mesh of + 1 
points. One defines the discrete approximation of i/M (■,£), z = 1 ? - - - , 7V C as 



u«(i)=[ M «(t),..., M « ) (i) 



(39) 



The space derivatives of u^) (t) are approximated as 

uW(t) w (pW) _1 QWuW(t), (40) 

where (pW, Q^), z = 1, • • • , iV c have similar properties as in the previous section. Diagonal matrices A^, 
GW(t) are defined as in the previous section. Similarly as in the previous section, the dependency of u and 
its derivatives on time is omitted for clarity. The semi-discretized equations are all of the same form as (30). 
The interface conditions (19)-(21) are added as penalty terms, resulting in 



N. 



(P«) 'qW ((a'^uW) - ^L A «G«(t)u« 



E 4:K p(<) )" 1 ( D " ) ) T K ) -^' ) )4 i) 



(41) 



-S(E(«o J) ) 2 ^-o 



(pw)" 1 ^) 



(i) 



where denotes the penalty term for the i-th cable associated with the boundary x = L^> : 

p«=a«(p«)- 1 ( u « (i) -0)e«, i = l,..., Ni -l 
for the sealed boundary conditions u^ N(i) = 0, and 

p W0 = ^c) JpW)" 1 ( D W)) T (u W0 _ 0)e ^) 

for the voltage clamp condition u^ N °^ = 0. These penalty terms can be all recast in the same form 



Pl 



(i) 



(pw)" 1 



7 yW pW 7 — 1 . . . AT- 



(42) 



(43) 



(44) 



Multiplying by (u^) P^ and using analog identities to the one derived in the previous section leads to 

2 



= -/i 

■ E 



aw u w 

E ff 8 u io («o 



/ GW(t)AW 

Cm, 



) _ 0) 
"o 



7V C 



(45) 



+ (^+M(ay w )>S rW 



W ^ .,(0 



I,-" ,JVc 
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Summing the iV c equalities (45) and choosing the penalty coefficients = —fx (o^ w ~^ , 

JO 



a$ = JL, i = l,... , N r 



and 

leads to the energy estimate 
2 



<iwr = ^-w ) 



1 Nc M 2 ^ M 2 ^ 

1^ k/XWu« =-^ A«u« -e 



lGW(t)A« 



u 



(0 



Defining the weighted norm 



|u||J = 2Ka(0u 



(i) 



2 = 1 



the following energy equality holds as well as the proposition 



Ml) + 2^ 



Au 



-eg,, 



0. 



(46) 
(47) 

(48) 

(49) 

(50) 



Proposition 5. The SBP-SAT scheme for solving the semi-discrete problem (41) associated with the initial 
boundary value problem (12), (19), (20), (21) is stable. 

Remark. Note that this energy estimate is the exact semi-discrete analog to the energy equality (26) derived 
m 

Section [2X2 



4. Order of convergence 

The method of manufactured solutions [21 is used to determine the discretization order. For the two 
problems mentioned above, a closed form solution is chosen and appropriate source terms are added to the 
initial boundary value problem of interest so that the closed form solution satisfies it. Time integration is 
done using the fourth-order explicit Runge-Kutta scheme and a time step of At = 10 -9 s, ensuring that the 
time discretization errors are negligible. The convergence rate is calculated as 




4-1. Cables with junction 

Three cables with junction at = 0, i = l, --3 are considered. Their respective cable radii are 
constant for each cable. The exact manufactured solutions u( % \x,t), i = 1, • • • ,3 and their corresponding 
initial boundary value problem are specified in Appendix C. The physical constants associated with the 
problems are reported in Tables [l]and|2l The exact solutions u^\x^t) at times t = and t = T are plotted 
in Figure [3] 

The solutions u^\x^\t) at t = and t = T are reported in Figure [3] 

The convergence results are reported in Tableland Figure [4j The slopes corresponding to each operators 
order are also shown in Figure [4] Second-order convergence is observed with second-order operators, third- 
order convergence with third-order operators, fourth-order convergence with fourth-order operators and 
fifth-order convergence with fifth-order operators. 

The CPU times associated with each computation are reported in Table [4j The CPU times corresponding 
to a relative error of the order of 10 -6 are also highlighted in bold font. For the same accuracy, the CPU time 
associated with the second-order operators is 2.5 times larger than the ones associated with the higher-order 
operators. 
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Specific membrane capacitance 


Cm 


10" 2 F.m~ 2 


Axial Resistivity 


Ri 


0.354 ohm.m 


Conductance 


9i 


1200 ohm-^rn- 2 


Conductance 


92 


360 ohm _1 .m -2 


Conductance 


9s 


3 ohm _1 .m -2 


Equilibrium potential 


E x 


0.115 V 


Equilibrium potential 


E 2 


-0.012 V 


Equilibrium potential 


E 3 


0.010613 V 



Table 1 Physical constants for the problems of interest 



Cable length 


lw 


0.05 m 


Cable length 


L( 2 ) 


0.05 m 


Cable length 


L( 3 ) 


0.063 m 


Cable radius 




0.476 x 10" 3 m 


Cable radius 




0.476 x 10" 3 m 


Cable radius 


aW(x) 


0.7556 x 10" 3 m 



Table 2 Geometric parameters associated with the junction problem 




Figure 3 Exact solutions for the junction problem 
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N 


order 2 


order 3 


order 4 


order 5 


32 


0.707003 


1.359596 


2.381819 


3.028414 


64 


0.969357 


2.641288 


3.757290 


5.092958 


128 


2.135884 


3.319466 


3.876460 


4.843739 


256 


2.099840 


3.097754 


4.138360 


4.951223 


512 


2.051456 


3.043053 


4.327195 


5.108929 



Table 3 Numerical convergence behavior for each operator choice 



Junction 




Mesh size 

Figure 4 Errors in function of mesh size and operators order for the junction problem 



N 


order 2 


order 3 


order 4 


order 5 


16 


33.05 


34.30 


34.26 


34.99 


32 


33.01 


34.54 


35.71 


35.74 


64 


34.88 


35.84 


38.97 


78.98 


128 


35.78 


39.48 


40.22 


86.54 


256 


39.47 


45.09 


47.87 


111.73 


512 


93.04 


125.99 


126.18 


190.01 



Table 4 CPU times (in seconds) for each operator choice in the junction problem 
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-1 

0.01 0.02 0.03 0.04 0.05 

x (m) 

Figure 5 Exact solutions for the cable with soma problem 



4.2. Cable and soma 

The exact manufactured solution u{x,t) and the corresponding initial boundary value problem are spec- 
ified in Appendix D. The physical constants associated with the problem are reported in Tables [T] and [5] 
The exact solutions u(x,t) at times t = and t = T are depicted in Figure [5] 



Cable length 


L 


0.05 m 


Cable radius 


a(x) 


0.476 x 10- 3 m 


Soma radius 


^soma 


2 x 10" 3 m 



Table 5 Geometric parameters associated with the cable with soma problem 

The convergence results are reported in Table [6] and Figure [6j in which slopes corresponding to the 
operators orders are also depicted. Second-order convergence can be observed using second-order operators, 
fourth order convergence is observed when using third- and fourth-order operators and fifth-order convergence 
for fifth-order operators. For this case, no real advantage can be seen in terms of using fourth-order operators 
over third-order operators. 

The CPU times associated with each computation are reported in Table (7) In that table, the CPU times 
corresponding to a relative error of the order of 5 x 10 _T are highlighted in bold font. Note that for the 
same accuracy, the CPU time associated with the second-order operators is 20 times larger than the ones 
associated with the high-order operators. 
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N 


order 2 


order 3 


order 4 


order 5 


32 


0.766683 


2.862748 


2.567863 


4.554248 


64 


1.345940 


3.386744 


3.004467 


4.927776 


128 


1.978948 


3.752389 


3.588062 


5.442091 


256 


2.005454 


3.962842 


4.090079 


5.691623 


512 


2.000201 


3.983892 


4.314146 


5.161448 



Table 6 Numerical convergence behavior for each operator choice 




Mesh size 

Figure 6 Errors in function of mesh size and operators order for the cable attached to soma problem 



N 


order 2 


order 3 


order 4 


order 5 


16 


17.03 


16.64 


16.03 


16.27 


32 


16.81 


16.78 


16.50 


16.67 


64 


20.25 


21.49 


20.47 


36.76 


128 


24.05 


24.80 


29.57 


42.30 


256 


39.76 


44.88 


49.14 


117.77 


512 


315.09 


294.97 


318.83 


551.21 



Table 7 CPU times (in seconds) for each operator choice in the cable attached to soma problem 
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5. Numerical applications 



5.1. Time discretization 

After space discretization by the approach developed in this paper, the Hodgkin-Huxley equations become 
a set of coupled ODEs of the form 

— Ai(w)u + bi(w,t) 

— A 2 (u)w + b 2 (u) 

where the discrete variables for the potential and gating variables are u G IR Ar+1 and w = [m, h, n] T G 
R 3(iv+i\ respectively. 

This set of coupled ODEs is then discretized in time, using the second-order staggered difference scheme 
initially developed by Hines [3]. The scheme is 

I - ^A 2 (u n )^ w n+ ^ = ^1 + ^A 2 (u n )^) w n ~i + Atb 2 (u n ) 



I 

2 



2 

At 



Ai(w n+ i)^ u n+1 = ^1+ ^Ax(w n+ ^ u n + Atbi(w n+ i,£ n+i )- 



Here u n « u(t n ) = u(nAt) and w n+ ^ « w(t n+ ^) = w + At 

5.#. Potential propagation in a cable attached to the soma 

The physical constants associated with this problem are provided in Table [TJ These correspond to the 
giant squid's axon studied by Hodgkin and Huxley in [1 j. The expressions for the gating functions are 
provided in Appendix A. 

Current with Gaussian distribution and maximum intensity I = 10 nA is injected at the top extremity 
of the cable depicted in Figure [2^a) and the cable potential time history for t = [0, 10] ms is recorded at the 
two extremities of the cable. The corresponding time histories obtained by applying the procedure developed 
in this paper with N — 32 and 5-th order differential operators and a time step size of dt = 10 -2 ms are 
reported in Figure [7[ The typical action potential characteristic shape [17] can be observed. 

Next, the action potential peak times obtained at the two extremities of the cable of interest are recorded 
and shown for various choices of order for the differential operators and various numbers of mesh points. 
The peak time at the injection points t p , depicted in Figure [7| are almost identical as it can be observed 
in Figure [8 However, for coarse meshes, second-order operators result in an inaccurate peak time t s p at the 
Soma, as observed in Figure [9j This results in inaccurate propagation times for these coarse meshes when 
second-order operators are used, whereas the propagations times are accurate with high-order operators even 



when using those coarse meshes, as shown in Figure 10 



5.3. Potential propagation in a dendritic network 

In order to illustrate the capability of the proposed approach to predict potential propagation in large 
dendritic networks, a tree with 15 branches attached to the soma, as depicted in Figure [IT] is considered. 
The dimensions of each branch in the tree are chosen according to Rail's 3/2 power law [2] and correspond 
to Rallpack's second configuration [22 . The cable lengths and radii at each level of the tree are reported in 
Table [8] All the other physical constants are identical as in the previous numerical experiments. Synaptic 
current input is modeled by a current injection of amplitude / = 500 nA for a duration of 1 ms at the 
extremity of each of the 8 extremal branches. Each branch is discretized using 31 nodes, resulting in 1860 
degrees of freedom. Fifth order difference operators are here used. The potential propagation is subsequently 
simulated in the entire tree in the time interval for t G [0, 0.2] s using a time step At = 10 -4 s. The entire 
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Figure 7 Action potential at the current injection point and the soma obtained using fifth-order operators 
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Figure 8 Action potential peak times f p at the current injection point 
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Figure 11 Dendritic tree connected to the soma 



simulation takes about 14 s CPU time when performed with MATLAB running on a Mac Book Pro 2.9 GHz 
Intel Core i7, 8GB 1600 MHz DDR3. 



Two current input configurations are considered. The first one is depicted in Figure 12 (a). In that figure, 
the interval of duration 1 ms, at which a given branch input is active, is represented by a line segment. The 
corresponding potential recorded at the soma is reported in Figure 12 (b). There are eight characteristic 
action potential spikes, each one corresponding to a specific extremal branch input. 

The second numerical experiment is carried out with the current input activity reported in Figure [l3| (a). 
In this case, the current input is limited to the time interval t G [0,0.1] s. The associated potential at the 
soma is presented in Figure [13] (b). Only four characteristic action potential spikes are propagated from the 
extremity branches of the tree to the soma. This numerical experimental emphasizes the nonlinear nature 
of the Hodgkin-Huxley equations. It also illustrates the concept of refractory period for action potentials 
observed in neurons [23]: after an action potential has been generated, the membrane is hyper-polarized 
and there is a minimal period of time after which another action potential can be created. In the second 
numerical experiment, although eight current inputs were injected, four failed to create full-amplitude spikes 
since their respective injection times were within the refractory periods of the four action potentials. A small 
spike can also be observed around t = 0.017 s, but it is not a full action potential. In the first experiment, 
however, the time between two current injections is longer, leading to eight spikes. 



Level 


Number of branches 


Branch length (m) 


Branch radius (m) 


1 


1 


32.0 x 10~ 6 


8.0 x 10~ 6 


2 


2 


25.4 x 10~ 6 


5.04 x 10~ 6 


3 


4 


20.16 x 10~ 6 


3.18 x 10~ 6 


4 


8 


16.0 x 10~ 6 


2.0 x 10~ 6 



Table 8 Dimensions of each branch in the dendritic tree 



17 



0.07 



8 
7 

CD 

CO o 

a 6 

^ 5 
m 

CD 

4 

o 

2 
1 



0.05 



0.1 

Time (s) 



0.15 



0.2 





0.06 


> 


0.05 








0.04 






CD 




o 


0.03 


Oh 




ma 


0.02 


o 










0.01 





-0.01 



V 



V 



V 



0.05 



0.1 

Time (s) 



0.15 



0.2 



(a) Current input activity for each of the extremal branches 
(each line segment is of length 1 ms) 



(b) Action potential at the soma 
Figure 12 Potential propagation in the network for the first current input configuration 



8 

7 

CO o 

a 6 

^5 

CD 

.% 4 

CD 

2 
1 



0.05 



0.1 

Time (s) 



0.15 



0.2 





0.07 




0.06 


> 


0.05 








0.04 






CD 




O 


0.03 






ma 


0.02 


o 




CO 






0.01 





-0.01 



V 



0.05 



0.1 0.15 

Time (s) 



(a) Current input activity for each of the extremal branches 
(each line segment is of length 1 ms) 



(b) Action potential at the soma 
Figure 13 Potential propagation in the network for the second current input configuration 
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6. Conclusions 

Stable and accurate high-order schemes based on Summation-By-Parts (SBP) form with the weak Si- 
multaneous Approximation Term (SAT) procedure for boundary and interface conditions are developed for 
the Hodgkin-Huxley equations. These schemes are shown in this paper to be easily applicable to neuronal 
systems constituted of multiple branches such as in the case of dendritic trees. Well-posedness of the con- 
tinuous problem as well as energy estimates are established for dendritic trees connected to the soma. The 
SBP-SAT procedure leads to an energy estimate and proof of stability of the discrete scheme. 

Convergence and order of accuracy are illustrated on model problems. Important gains in CPU times 
for a given level of accuracy are also demonstrated when high-order operators are used. The superiority of 
high-order operators is also illustrated on a simple problem constituted of a cable connected to the soma. 
Finally, the proposed approach is applied to a realistic dendritic tree configuration constituted of 15 cables 
connected to the soma. On that configuration, the nonlinearity properties of the Hodgkin-Huxley equations 
are numerically demonstrated. 

Appendix A: expressions for the gating functions 

The expressions for the coefficients o^ m7 

/? m , f3 h and f3 n were determined for the giant squid axon 

by Hodgkin and Huxley pQ as: 

, v 10 5 (0.025 -ii) r> / \ A ^_ 

e o.oi — 1 

10 3 

a h (u) = 70e~^, P h {u) . 03 _, (52) 

e o.oi -f 1 

, >. 10 4 (0.01 -U) r> / \ 

<*n(u) = — -, Pn{u) = 125e o.os . 

e o.oio — 1 



Appendix B: proof of proposition 1 

The differential equations associated with m{x,t) are a set of uncoupled ODEs for each x. As a result, 
it is sufficient to consider one of those ODE only for xo G [0, L]: 

m t (xo,t) = a m (u(xo,t))(l - m(x ,t)) - /3 m (u(xo,t))m(x ,t), t G [0,T]. (53) 

To simplify the notations, in the rest of this section, the dependency on x is dropped: a m (u(xo,t)) = 
a m (u(i)), f3 m (u(x ,t)) = P m (u(t)) and m(x ,t) = m(t). 

From the expressions of a m and /3 m given in Appendix 1, it can be shown that 

a m (0)>0, /3 m (l)>0. (54) 

These properties also hold for o^, a n: fih and f3 n . 

An initial condition m(0) G [0, 1] is then considered. Assume that there exists t\ G (0, T] such that 
m(ti) > 1. As ra(-) is continuous, this means that there exists < to < t\ < T such that m(to) = 1 and 
m(to + e) > 1, VO < e < t\ - t . Since 

dim 

— (t ) = a m (u(to))(l - 1) - /3 m (u(t ))l = -Pm(u(t )) < 0, (55) 

this is a contradiction. Hence m(t) < 1, \/t G [0,T]. 

Similarly, assume now that there exists t\ such that m(t\) < 0. As m(t) is continuous, that means that 
there exists < to < t\ such that m(to) = and m(to + e) < 0, VO < e < t\ — to. However, 

^(t ) = a m (u(to))(l - 0) - p m (u(t ))0 = a m (u(t )) > 

which leads to a contradiction as well. 

In conclusion < m(t) < 1, Vt G [0, T]. This property also holds for h{t) and n(i), respectively. 
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Appendix C: manufactured solution for the junction problem 



A case with 3 branches having a junction at #W = o, i = 1, • • • , 3 is considered. The geometric parameters 
of the problem are a (1) (x (1) ) = a (2) (x (2) ) = 2" 2/3 a (3) (x (3) ) = a , where e [0,L^], and their lengths are 
given by = L, i = 1,2 and L^ 3 ^ = 2 1 / 3 L. The following manufactured solutions are considered, for a 
problem without current injection, that is I(t) = 0, t G [0,T], where T = 10 -5 s: 



*W(a;W,t) =exp 



^ 3 )(^ 3 \t) = exp 



m 



i«(x«,t) 



The solutions ( 56 ) satisfy the equations 



exp 




( 1 3 


exp 


;H 


( i 3 


1, i 


= 1, 


••• ,3, 


1, i 


= 1, 


••• ,3, 


1, i 


= 1, 


••• ,3, 


], * = 


= 1,- 


• ,3. 




sin 



sin 



c 



37nr« 
V 2L 

2 4 /3L j 



1,2 



(56) 



m[ i} = « m (ttW)(l - mW) - fc(« (,) )m (,) + 



(i) 



a h («W)(l - - /3 fc (« W )/i (i) + if } ( U «) 
<*„(«<*>)(! - ra«) - ^„(«W)nW + („(*)), 



(57) 



where 



1 3 



F«(i) = /?«(««) 

^°(t) = /#V°) 
F«(t) = /3«( U «) 



(58) 



and the boundary conditions are 



u (i) (0,t) 


= 0, 




= 0, 




= 0, i 




= exp 



(59) 



3tt 
2L 



The numerical values chosen for each physical constant are the same as the ones reported in Table [T] and [2] 
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Appendix D: manufactured solution for the cable with soma problem 



A problem with constant cable radius a(x) = ao, x G [0,L] and without current injection, that is 
I(t) = 0, t G [0,T], where T = 1(T 5 s: 





= exp 


m(x, t) 


= 1 


ft(x, t) 


= 1 


n(x, t) 


= 1 



cos 



/3x 



(60) 



where (x,t) G [0,L] x [0,T] and /3 satisfies ^ 
The solutions ([60]) satisfy the equations 



77a L ' 



/i 



1 



(a(x) 2 u x ) x - —g(m, ft, n)u + /(ra, ft, n, x, t) + 



a(x) v v y C n 
m t = a m (u)(l - m) - P m (u)m + F m (w) 
ft £ = a,»(l - ft) - p h {u)h + i^) 
n* a n (w)(l - n) - /3 n (u)n + F n (u) 



(61) 



where 



F u {t) 



C, 



1 3 

— Hft^i 



i=i 



= 0m («) 

F„(i) = /3 n (u) 



(62) 



and the boundary conditions 
u*(0) = 

u t (L) = -r]a(L) 2 u x (L) 

where 



1 #(m(L), ft(L), n(L)) + /(m(L), ft(L), n(L), L, *) + F 6 (t) iM 



l 3 



(64) 



The numerical values chosen for each physical constant are provided in Tables [T] and [5] 
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